Numerical analysis of the master equation 



Ronald Dickman 1 

Departamento de Fisica, ICEx, Universidade Federal de Minas Gerais, Caixa Postal 702, 

30161-970, Belo Horizonte - MG, Brasil 

(February 1, 2008) 

Abstract 



Applied to the master equation, the usual numerical integration meth- 
ods, such as Runge-Kutta, become inefficient when the rates associated 
with various transitions differ by several orders of magnitude. We intro- 
duce an integration scheme that remains stable with much larger time 
increments than can be used in standard methods. When only the sta- 
tionary distribution is required, a direct iteration method is even more 
rapid; this method may be extended to construct the quasi- stationary 
distribution of a process with an absorbing state. Applications to birth- 
and-death processes reveal gains in efficiency of two or more orders of 
magnitude. 
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The master equation is the basic tool for describing Markovian stochastic pro- 
cesses on a discrete state space, in continuous time Despite its central role 
in stochastic analysis in the physical and biological sciences, only a handful of ex- 
actly soluble examples are known. Thus the need to integrate the master equation 
numerically arises frequently, and with it the issue of computational efficiency 

Consider a Markov process with transition rates W n m from state m to state n. 
In the master equation (ME), 

Pn = W n,n'Pn' ~ Pn Y W n',n (n = 0, N), (1) 
n' n' 

the factor J2 n ' ^n'.n = w n multiplying p n can become large in a typical birth-and- 
death process, where W n > >n ocnor some higher power of n. This poses a problem for 
numerical integration via the usual discretization schemes, such as the Runge-Kutta 
method (RKM). Since instabilities appear when |w n At| > 1, we must use a small 
time increment, At ~ l/(max n w n ), and convergence is slow. With some (but not 
all) of the transition rates large, the ME is, in effect, a stiff system of differential 
equations, requiring special numerical treatment [[| . In this note I introduce numer- 
ical integration schemes for the ME that are efficient when the transition rates vary 
over a wide range, and an iteration method that eliminates the need for step-by-step 
integration, when only the stationary (or quasi-stationary) distribution is required. 
Simple but detailed examples are used to illustrate the methods; further applica- 
tions will be reported elsewhere |3J]. (We consider stationary Markov processes, i.e., 
time-independent rates, although the method is not limited to this class of problem.) 

We begin by writing the ME in the form 

Pn = -W n p n + T n , (2) 

where 

r n (t) =J2 W n,n'Pn'(t) . (3) 
n' 

Note that p n does not appear in the sum for r n , since W n ,n = 0. Both w n and r n 
are nonnegative. (In fact w n is zero only if state n is absorbing.) Integrating Eq. 
(0) we have 

Pn(t) = e-^Vn(O) + f dt'e-^-^r^t') . (4) 

Jo 

This is only a formal solution, since we need the p n (t) to evaluate r ra , but it is a 
useful starting point for approximations. If we adopt a time increment At such that 
r n (t) is approximately constant over this interval, then we have 
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Pn(At) 



-w n At 



Pn(0) + 



-w n At 



r n (0) 



(5) 



This relation can then be iterated: we use p n (At) to evaluate r n at the start of the 
second interval, and thereby find p„(2At), and so on. This simple formal integration 
(FI) scheme, analogous to Euler's method for direct numerical integration, already 
represents a significant advantage over the usual approaches when some of the w n are 
large. The reason is that the exponential factor is already included in the solution, 
whereas in the usual discretizion methods it has to be constructed term by term, in 
powers of w n At. 

Suppose that the ME of interest possesses a unique, stable stationary solution p n . 
For a stationary solution, p n = r n /w n , where r n is given by Eq. (||) with p n = p n . 
Setting p n = p n on the r.h.s. of Eq. (Q) we immediately see that p n is a stationary 
solution of our iteration method. To investigate the linear stability of the stationary 
solution, let p n = p n + 5 n . According to the ME, 

S n = W n,n'8 n > ~ S n ^ W n > ,„ , (6) 
n' n' 



or, in matrix notation, 



±S = 8S (7) 



where S nim = W n ^ m for n ^ m and S n , n = —w n . Since, by hypothesis, the stationary 
solution 5 = is stable, S has one zero eigenvalue and all others negative. 

If we insert p n = p n + 5 n in the r.h.s. of Eq. @ (so that r n = r n + J2 n ' W ntn '5 n >), 
we find 

S(t + At) = T5(t) (8) 



where T„ jm = (1 - e WnAt )W nim /w n for n ^ m, and T ni „ = e WnAt . Expanding to 
first order in At, we see that Eq. (^) applies in this well, so the stationary 

solution is stable if it is so for the original ME. 

We refer to the integration method embodied in relation Eq. (|5]) as a "first-order 
FI scheme" since the error in the solution p n (t) is proportional (for fixed t) to At. 
[In each step, the error incurred in treating r n as constant is O(At) 2 (since in fact 
r n (At) = r n (0) + r' n (0)At + •••), and we require N = t/At steps.] An obvious way 
to improve accuracy, analogous to going from Euler's method to the midpoint (or 
second-order Runge-Kutta) method, is to replace the assumption of a constant r n 
on the interval [t, t + At] with a linear approximation. To do this we first use Eq. 
(|) to find p n {t + At), and then estimate the r n (t + At) using these values. Then 
we form the linear approximation 
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r n ( s )~r n (t) + (s-t)r' n , (t < s < t + At) , (9) 



where 

, r n (t + At) - r n (t) 



(10) 

At v ; 



Inserting this in the formal solution, Eq. (HI), we find 



Pn {t + At) ~ e- w - At p n {t) + 



1 - e 



r„(t) 



+ 



UK 



At 



1-e 



-w„At~ 



— (11) 



(Note that we do not use the expression for dr n /dt that follows from the master 
equation in place of r' n , since it does not, in general, represent the behavior of r n 
over a substantial time interval.) The total error associated with the second-order 
FI procedure is oc (At) 2 . Further improvement is clearly possible, by introducing 
higher order polynomial approximations to the r n , but the second-order scheme is 
quite adequate for the applications considered here. The error estimates apply for 
a fixed, finite time; the error for t — > oo is zero, if the system possesses a unique 
stationary distribution. 

A tremendous simplification and speedup is possible if one only requires the 
stationary probability distribution. The relation p n = r n /w n suggests an iterative 
procedure of the form 

r 

p' n = ap n + (1 - a)— (12) 



where < a < 1 is a parameter. We expect p' n to converge to Ap n , where the factor 
A depends on a and on the initial distribution. 

Some insight into the choice of a is afforded by the simple example of a two- 
state system with transition rates W\ t o = A and Wq,i = The evolution of the 
probability distribution is given (in matrix notation) by p' = Mp, where 



M 



( 



V(l 



:i-«)f 



a)* 



(13) 



Matrix M has eigenvalues 1 and 2a — 1. Thus the iterative procedure converges to 
a multiple of the stationary distribution, (//, A) for < a < 1, and is instantaneous 
for a = 1/2. 

A similar analysis of the three-state processes (n = 0, 1, 2), with rates 
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a) W n+ltn = 7, W n - 1>n = 1, 



and 



b) W n+hn = 7, W 2 , = 7 2 , and W n ^ m = 1 for n < m, 

yields a = 1/3 as the optimal choice. (By optimal choice we mean the value that 
minimizes maxj |a>j|, the uii being the eigenvalues of M, excluding, of course Uq = 1, 
associated with the stationary distribution.) These results suggest that for processes 
with a large number of states, even smaller values of a may be advantageous. The 
numerical examples discussed below support this notion. 

A simple modification of the iterative scheme, Eq. fll2|) , generates the quasi- 
stationary distribution (if such exists) of a Markov processs with an absorbing state. 
Consider a process on the states n — 0, 1, 2, ...N with n = absorbing, i.e., W n ,o = 
for all n, while Wo t7l for at least one n > 0. (It seems reasonable to exclude manifestly 
transient states, in other words, we assume that for each n, W n>m > for at least 
one m.) In this case the stationary state is p n = 8 nj o, but it is possible that the 
probability distribution, conditioned on survival, attains a stationary form, that is, 
for long times p n (t) — > C{t)q n (n > 0), where the q n are time- independent. Such 
quasi- stationary distributions arise in birth-and-death processes with saturation, for 
example the Malthus-Verhulst process or the contact process [|J. 

The defining feature of the quasi-stationary distribution is that relative proba- 
bilities Pn(t)/p m (t) (n, m > 0) are constant, or equivalently, p n /Pn = k, constant 
and independent of n for n > 1. Suppose the probability distribution has attained 
a quasi-stationary form at some time t, and that at this instant the distribution 
is normalized so: J2n>iPn(t) = 1. Using Eq. @ we have r n = (k + w n )p n , and 
summing on n > 1 yields 

K = Pn = ~ J2 W 0,nPn = -Tq ■ (14) 
n>l n>\ 

(ro is the decay rate of the survival probability in the quasi-stationary regime.) Thus 
we find that in the quasi- stationary state, 

Pn = Tn , {n > 1). (15) 
w n - r 

This relation suggests that we iterate as follows: 

T 

Pn = a Pn + {l-a) • (16) 



The new distribution p' should be normalized after each iteration, since Eq. QT5j ) 
assumes this property. [If r = the process of course possesses a true station- 
ary distribution and Eq. ( |l6l) reduces to Eq. flO).] We have verified that this 
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scheme converges to the quasi-stationary distribution much more rapidly than via 
integration of the ME. 

As a first example we consider the nonlinear one-step process (on states n = 
1, 2, 3, ...), with nonzero transition rates: 

Wn-i,n =n(n-l) , W n+ltn = X (17) 



This represents the coagulation process A + A — > A with addition of species A at 
rate A, in a well-stirred system. In Fig. 1 we compare the mean population size 
(n)t as furnished by integration of the ME via a fourth-order RKM (solid line) with 
the second-order FI scheme, Eq. (PI), (open circles). [Here A = 100, p n (0) = 5„ i, 
and we have set p n (t) = for n > 100.] The results are nearly identical, but the 
RKM requires and integration step of At < 2 x 10~ 4 for stability, whereas the FI 
integration uses At = 0.002. The latter yields (n) t with a relative error of less than 
0.8%, the largest deviation occurring at relatively short times. (As noted above, the 
stationary values are identical in all cases.) To maintain fidelity to the solution at 
short, as well as long times, one may use the FI with an adjustable time step (+ 
in Fig. 1), such that r'*At is small, where r'* = max n {r' n } . For example, using the 
second-order FI with a time step of At = 0.0005 + .02/(1 + r'*), the relative error is 
reduced to < 0.06%, while the average time step over the region of interest is 0.0024. 

Fig. 2 shows the stationary probability distribution for A = 400. The solid line 
represents the RKM result, while the points come from iteration of Eq. ([12]) . The 
distributions are identical to within one part in 10 5 . In this case, the distribution 
converges (such that the error in (n) t is < 10~ 10 ), after about 900 iterations, when 
we use a = 1/2. As a is reduced, the number of iterations required falls steadily; in 
fact the procedure converges even for a = 0, after only 460 steps. 

Our second example is a multi-step generalization of the Malthus-Verhulst pro- 
cess, with transition rates: 

f n[l+/i(n-l)]e 7(m " n+1) , m<n 
W m , n = < Xne^ n - m+1 ^ , m>n (18) 

[ , m = n 



for n, m = 0,1,2,.... (Thus transition rates between states n and m fall off oc 
exp[— j\m — n|]. We choose this example to show that the methods proposed here 
are not limited to one-step processes, which admit a relatively simple analysis 
In this case n = is absorbing. Figure 3 compares the quasi-stationary distribution 
for A = 3, /i = 0.1, and 7 = 1, obtained via the RKM (with At = 10~ 4 ; larger 
values produce instability), with iteration of Eq. (JT6]). As is clear from the plot, 
the results are identical in every detail. The RKM requires about 1270 seconds of 
cpu on a DEC-alpha workstation to converge to the quasi-stationary distribution; 
the FI method requires about 50 sec. Iteration of Eq. (|T6|), (with a = 0.5), by 
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contrast, converges in 9 seconds (about 1600 steps). Reducing parameter a to 0.1 
yields convergence in just 900 steps (5 sec. cpu), and in 815 steps for a = 0. 



The steady decrease in computation time as we reduce a leads one to ask whether 
it is possible to use a negative. Iteration of Eq. (|T6| ) with a < does not work, as 
negative probabilities are generated, but we can exclude these by writing 



Pn 



max 



0, ap n + (1-a)- 



W r . 



■r 



(19) 



This functions even for negative a, converging, in the present case, after about 630 
steps when a = —0.3. For a < —0.4, however, the scheme does not converge. On 
the other hand, in the first example [i.e., rates given by Eq. flT7D], using a < does 
not offer any advantage over a = 0. 

We conclude from these and other examples with large numbers of states [[|, 
that a ~ generally yields rapid convergence, but that some amount of experiment- 
ing (for example, with negative a values), may reduce computation time further, 
in particular cases. But even without such optimization, the iterative schemes of 
Eqs. (0) and (|T6"|) yield economies in computation time of two or more orders of 
magnitude, compared with conventional integration methods. As the number of 
variables and/or the disparity in the magnitudes of the transition rates increases, 
direct iteration and the FI scheme become ever more valuable. 
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FIGURE CAPTIONS 



FIG. 1. Mean population (n) versus time in the coagulation process, Eq. (|T7|), with 
A = 100, using Runge-Kutta integration (solid line), the second-order FI scheme 
with a fixed time step (circles), and with a variable time step, as described in the 
text (+). 



FIG. 2. Stationary probability distribution in the coagulation process with A = 400. 
Solid line: RKS; +: direct iteration of Eq. (|T2"D. 



FIG. 3. 

Verhulst process 
+: direct iteration of Eq. (|l( 



Quasi- stationary probability distribution in the generalized Malthus- 

0.1, and 7 



Eq. (|T8]), with A = 3, /i 



1. Solid line: RKS; 
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